############# Blame Experiment, Study 1, Lucid and TESS samples ##############
# Author: Ryan Kennedy

# Load needed packages
library(tidyverse)
library(here)
library(arm)
library(tidybayes)
library(effsize)
library(skimr)
library(broom)
library(patchwork)
library(stargazer)
library(ggpubr)

# Load combined data
tess2 <- read_csv(here("Data","combined_study_1.csv"), col_types = cols(tia = col_double()))
tess2 <- tess2 %>% mutate(treatment = factor(treatment, levels = c("j", "a_a", "a_d",
                                                                   "h_a", "h_d",
                                                                   "ha_a", "ha_d")))

# Create bar chart chart of blame allocation in different treatments
descplotc <- tess2 %>%
  gather(allocate, amount, jblame:ablame) %>%
  group_by(treatment, allocate) %>%
  summarize(mamount = mean(amount, na.rm = T)) %>%
  ungroup() %>%
  mutate(allocate2 = case_when(allocate == "ablame" ~ "Algorithm",
                               allocate == "hblame" ~ "Human",
                               allocate == "jblame" ~ "Judge"),
         treatment2 = case_when(treatment == "j" ~ "Control",
                                treatment == "a_a" ~ "Agrees with algorithm",
                                treatment == "a_d" ~ "Disagrees with algorithm",
                                treatment == "h_a" ~ "Agrees with human",
                                treatment == "h_d" ~ "Disagrees with human",
                                treatment == "ha_a" ~ "Agrees with algorithm and human",
                                treatment == "ha_d" ~ "Disagrees with algorithm and human")) %>%
  mutate(treatment2 = factor(treatment2, levels = c("Control","Agrees with algorithm",
                                                    "Disagrees with algorithm",
                                                    "Agrees with human",
                                                    "Disagrees with human",
                                                    "Agrees with algorithm and human",
                                                    "Disagrees with algorithm and human")),
         allocate2 = factor(allocate2, levels = c("Judge","Algorithm", "Human"))) %>%
  ggplot() +
  geom_bar(aes(x = allocate2, y = mamount, fill=factor(ifelse(allocate2=="Judge","Highlighted","Normal"))), stat = "identity", position = "dodge") +
  facet_grid(~treatment2, labeller = label_wrap_gen(width = 2, multi_line = TRUE)) + 
  scale_fill_manual(name = "area", values=c("red4","grey50")) +
  labs(x = "Actor", y = "Average level of blame", title = "(A) Blame allocation") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), legend.position = "none")
descplotc

# Create boxplot of blame allocation in different treatments (Note: warning is because under some conditions there is not an option of algorithm and/or human blame by design)
descplotbox <- tess2 %>%
  gather(allocate, amount, jblame:ablame) %>%
  mutate(allocate2 = case_when(allocate == "ablame" ~ "Algorithm",
                               allocate == "hblame" ~ "Human",
                               allocate == "jblame" ~ "Judge"),
         treatment2 = case_when(treatment == "j" ~ "Control",
                                treatment == "a_a" ~ "Agrees with algorithm",
                                treatment == "a_d" ~ "Disagrees with algorithm",
                                treatment == "h_a" ~ "Agrees with human",
                                treatment == "h_d" ~ "Disagrees with human",
                                treatment == "ha_a" ~ "Agrees with algorithm and human",
                                treatment == "ha_d" ~ "Disagrees with algorithm and human")) %>%
  mutate(treatment2 = factor(treatment2, levels = c("Control","Agrees with algorithm",
                                                    "Disagrees with algorithm",
                                                    "Agrees with human",
                                                    "Disagrees with human",
                                                    "Agrees with algorithm and human",
                                                    "Disagrees with algorithm and human")),
         allocate2 = factor(allocate2, levels = c("Judge","Algorithm", "Human"))) %>%
  ggplot() +
  geom_boxplot(aes(x = allocate2, y = amount, fill=factor(ifelse(allocate2=="Judge","Highlighted","Normal"))), notch = T) +
  facet_grid(~treatment2, labeller = label_wrap_gen(width = 2, multi_line = TRUE)) + 
  scale_fill_manual(name = "area", values=c("red4","grey50")) +
  labs(x = "Actor", y = "Average level of blame", title = "(A) Blame allocation") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), legend.position = "none")
descplotbox

# Run basic model for effect of treatments on blame on judge
cm1 <- lm(jblame ~ treatment, data = tess2);summary(cm1)
glance(cm1)
tidy(cm1)

# Plot the results from the models, with coefficient estimates and 95% confidence intervals
modplotc <- data.frame(coef(sim(cm1, n.sims = 1000))) %>%
  dplyr::select(-X.Intercept.) %>%
  gather(treatment, sims, treatmenta_a:treatmentha_d) %>%
  group_by(treatment) %>%
  summarize(meffect = mean(sims),
            sdeffect = sd(sims)) %>%
  ungroup() %>%
  mutate(ehi = meffect + 2*sdeffect,
         elo = meffect - 2*sdeffect) %>%
  mutate(treatlabs = case_when(treatment == "treatmenta_a" ~ "Agrees with algorithm",
                               treatment == "treatmenta_d" ~ "Disagrees with algorithm",
                               treatment == "treatmenth_a" ~ "Agrees with human",
                               treatment == "treatmenth_d" ~ "Disagrees with human",
                               treatment == "treatmentha_a" ~ "Agrees with algorithm and human",
                               treatment == "treatmentha_d" ~ "Disagrees with algorithm and human")) %>%
  mutate(treatlabs = factor(treatlabs, levels = c("Disagrees with algorithm and human",
                                                  "Agrees with algorithm and human",
                                                  "Disagrees with human",
                                                  "Agrees with human",
                                                  "Disagrees with algorithm",
                                                  "Agrees with algorithm"))) %>%
  ggplot() +
  geom_errorbar(aes(x = treatlabs, ymin = elo, ymax = ehi, 
                    color=factor(ifelse(treatlabs=="Agrees with algorithm" | 
                                          treatlabs=="Disagrees with algorithm" |
                                          treatlabs=="Agrees with algorithm and human" |
                                          treatlabs=="Disagrees with algorithm and human","Highlighted","Normal")))) +
  geom_point(aes(x = treatlabs, y = meffect)) + 
  scale_color_manual(name = "area", values=c("royalblue3","grey60")) +
  geom_hline(yintercept = 0) + coord_flip() + ylim(-0.025, 0.15) +
  labs(x = "Treatment", y = "ATE", title = "(B) Average Treatment Effects (ATE)") + theme_bw() +
  theme(legend.position = "none")
modplotc

# Combine the plots
ggarrange(descplotbox, modplotc, nrow = 2)

# R-run analyses with for Lucid and TESS samples separately
cm2 <- lm(jblame ~ treatment, data = tess2[tess2$sample == "Lucid",]);summary(cm2)
cm3 <- lm(jblame ~ treatment, data = tess2[tess2$sample == "TESS",]);summary(cm3)

# Create table of main results
stargazer(cm1,
          dep.var.labels = c("Blame on judge"),
          covariate.labels = c("Agrees with algorithm",
                               "Disagrees with algorithm",
                               "Agrees with human",
                               "Disagrees with human",
                               "Agrees with algorithm and human",
                               "Disagrees with algorithm and human"),
          type = "html",
          out = here("Tables", "main_model.doc"))

# Create table of results for separate samples
stargazer(cm2, cm3, 
          dep.var.labels = c("Blame on judge (Lucid)", "Blame on judge (TESS)"),
          covariate.labels = c("Agrees with algorithm",
                               "Disagrees with algorithm",
                               "Agrees with human",
                               "Disagrees with human",
                               "Agrees with algorithm and human",
                               "Disagrees with algorithm and human"),
          type = "html",
          out = here("Tables", "split_blame_model.doc"))

# Run combined models for algorithm and human blame
cm_a <- lm(ablame ~ treatment, data = tess2);summary(cm_a)
cm_h <- lm(hblame ~ treatment, data = tess2);summary(cm_h)

# Create table of algorithm and human blame
stargazer(cm_a, cm_h, 
          dep.var.labels = c("Blame on algorithm", "Blame on human"),
          covariate.labels = c("Disagrees with algorithm",
                               "Disagrees with human",
                               "Agrees with algorithm and human",
                               "Disagrees with algorithm and human"),
          type = "html",
          out = here("Tables", "other_blame_model.doc"))

# Run combined model with covariates
cm1c <- lm(jblame ~ treatment + female + age + pid7 + incomec + white + educ5, data = tess2);summary(cm1c)

# Create table of results with covariates for SI
stargazer(cm1c, 
          dep.var.labels = c("Blame on judge"),
          covariate.labels = c("Agrees with algorithm",
                               "Disagrees with algorithm",
                               "Agrees with human",
                               "Disagrees with human",
                               "Agrees with algorithm and human",
                               "Disagrees with algorithm and human",
                               "Female", "Age", "Party ID", 
                               "Income", "White", "Education"),
          type = "html",
          out = here("Tables", "covariate_model.doc"))

# Plot these results
modplotcc <- data.frame(coef(sim(cm1c, n.sims = 1000))) %>%
  dplyr::select(treatmenta_a:treatmentha_d) %>%
  gather(treatment, sims, treatmenta_a:treatmentha_d) %>%
  group_by(treatment) %>%
  summarize(meffect = mean(sims),
            sdeffect = sd(sims)) %>%
  ungroup() %>%
  mutate(ehi = meffect + 2*sdeffect,
         elo = meffect - 2*sdeffect) %>%
  mutate(treatlabs = case_when(treatment == "treatmenta_a" ~ "Agrees with algorithm",
                               treatment == "treatmenta_d" ~ "Disagrees with algorithm",
                               treatment == "treatmenth_a" ~ "Agrees with human",
                               treatment == "treatmenth_d" ~ "Disagrees with human",
                               treatment == "treatmentha_a" ~ "Agrees with algorithm and human",
                               treatment == "treatmentha_d" ~ "Disagrees with algorithm and human")) %>%
  mutate(treatlabs = factor(treatlabs, levels = c("Disagrees with algorithm and human",
                                                  "Agrees with algorithm and human",
                                                  "Disagrees with human",
                                                  "Agrees with human",
                                                  "Disagrees with algorithm",
                                                  "Agrees with algorithm"))) %>%
  ggplot() +
  geom_errorbar(aes(x = treatlabs, ymin = elo, ymax = ehi, 
                    color=factor(ifelse(treatlabs=="Agrees with algorithm" | 
                                          treatlabs=="Disagrees with algorithm" |
                                          treatlabs=="Agrees with algorithm and human" |
                                          treatlabs=="Disagrees with algorithm and human","Highlighted","Normal")))) +
  geom_point(aes(x = treatlabs, y = meffect)) + 
  scale_color_manual(name = "area", values=c("royalblue3","grey60")) +
  geom_hline(yintercept = 0) + coord_flip() + ylim(-0.025, 0.15) +
  labs(x = "Treatment", y = "ATE", title = "Average Treatment Effects (ATE)") + theme_bw() +
  theme(legend.position = "none")
modplotcc

# Check some interactions with demographics to estimate CATEs
cm1cf <- lm(jblame ~ treatment + treatment * female + female + age + pid7 + incomec + white + educ5, data = tess2);summary(cm1cf)
cm1ca <- lm(jblame ~ treatment + treatment * age + female + age + pid7 + incomec + white + educ5, data = tess2);summary(cm1ca)
cm1cw <- lm(jblame ~ treatment + treatment * white + female + age + pid7 + incomec + white + educ5, data = tess2);summary(cm1cw)
cm1ce <- lm(jblame ~ treatment + treatment * educ5 + female + age + pid7 + incomec + white + educ5, data = tess2);summary(cm1ce)
cm1ci <- lm(jblame ~ treatment + treatment * pid7 + female + age + pid7 + incomec + white + educ5, data = tess2);summary(cm1ci)
cm1ct <- lm(jblame ~ treatment + treatment * tia + female + age + pid7 + incomec + white + educ5 + tia, data = tess2[tess2$sample == "Lucid",]);summary(cm1ct)

# Create table of results
stargazer(cm1cf, cm1ca, cm1cw, cm1ce, cm1ci, cm1ct, 
          dep.var.labels = c("Blame on judge"),
          covariate.labels = c("Agree with algorithm",
                               "Disagree with algorithm",
                               "Agree with human",
                               "Disagree with human",
                               "Agree with algorithm and human",
                               "Disagree with algorithm and human",
                               "Trust in automation",
                               "Female",
                               "Age",
                               "Ideology", "Income", "White",
                               "Agree with algorithm * education",
                               "Disagree with algorithm * education",
                               "Agree with human * education",
                               "Disagree with human * education",
                               "Agree with algorithm and human * education",
                               "Disagree with algorithm and human * education",
                               "Education",
                               "Agree with algorithm * female",
                               "Disagree with algorithm * female",
                               "Agree with human * female",
                               "Disagree with human * female",
                               "Agree with algorithm and human * female",
                               "Disagree with algorithm and human * female",
                               "Agree with algorithm * age",
                               "Disagree with algorithm * age",
                               "Agree with human * age",
                               "Disagree with human * age",
                               "Agree with algorithm and human * age",
                               "Disagree with algorithm and human * age",
                               "Agree with algorithm * white",
                               "Disagree with algorithm * white",
                               "Agree with human * white",
                               "Disagree with human * white",
                               "Agree with algorithm and human * white",
                               "Disagree with algorithm and human * white",
                               "Agree with algorithm * ideology",
                               "Disagree with algorithm * ideology",
                               "Agree with human * ideology",
                               "Disagree with human * ideology",
                               "Agree with algorithm and human * ideology",
                               "Disagree with algorithm and human * ideology",
                               "Agree with algorithm * TIA",
                               "Disagree with algorithm * TIA",
                               "Agree with human * TIA",
                               "Disagree with human * TIA",
                               "Agree with algorithm and human * TIA",
                               "Disagree with algorithm and human * TIA"),
          type = "html",
          out = here("Tables", "interaction_model.doc"))

#################
# Check for moderation based on covariates
#################

# Control -- only for TIA, which is not in the TESS survey
m_c_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
              data = tess2[tess2$treatment == "j",]);summary(m_c_tia)
# Agree with algorithm
m_a_a_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
                data = tess2[tess2$treatment == "a_a",]);summary(m_a_a_tia)
# Disagree with algorithm
m_a_d_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
                data = tess2[tess2$treatment == "a_d",]);summary(m_a_d_tia)
# Agree with human
m_h_a_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
                data = tess2[tess2$treatment == "h_a",]);summary(m_h_a_tia)
# Disagree with human
m_h_d_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
                data = tess2[tess2$treatment == "h_d",]);summary(m_h_d_tia)
# Agree with both
m_ha_a_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
                 data = tess2[tess2$treatment == "ha_a",]);summary(m_ha_a_tia)
# Disagree with both
m_ha_d_tia <- lm(jblame ~ tia + female + age + pid7 + incomec + white + educ5, 
                 data = tess2[tess2$treatment == "ha_d",]);summary(m_ha_d_tia)

# Calculate ATMEs for TIA
tia_a_a_coef <- coef(m_a_a_tia)[2] - coef(m_c_tia)[2]
tia_a_a_se <- sqrt(summary(m_a_a_tia)$coefficients[2,2]^2 + summary(m_c_tia)$coefficients[2,2]^2)

tia_a_d_coef <- coef(m_a_d_tia)[2] - coef(m_c_tia)[2]
tia_a_d_se <- sqrt(summary(m_a_d_tia)$coefficients[2,2]^2 + summary(m_c_tia)$coefficients[2,2]^2)

tia_h_a_coef <- coef(m_h_a_tia)[2] - coef(m_c_tia)[2]
tia_h_a_se <- sqrt(summary(m_h_a_tia)$coefficients[2,2]^2 + summary(m_c_tia)$coefficients[2,2]^2)

tia_h_d_coef <- coef(m_h_d_tia)[2] - coef(m_c_tia)[2]
tia_h_d_se <- sqrt(summary(m_h_d_tia)$coefficients[2,2]^2 + summary(m_c_tia)$coefficients[2,2]^2)

tia_ha_a_coef <- coef(m_ha_a_tia)[2] - coef(m_c_tia)[2]
tia_ha_a_se <- sqrt(summary(m_ha_a_tia)$coefficients[2,2]^2 + summary(m_c_tia)$coefficients[2,2]^2)

tia_ha_d_coef <- coef(m_ha_d_tia)[2] - coef(m_c_tia)[2]
tia_ha_d_se <- sqrt(summary(m_ha_d_tia)$coefficients[2,2]^2 + summary(m_c_tia)$coefficients[2,2]^2)

# Control -- for all the non-TIA variables
m_c <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
          data = tess2[tess2$treatment == "j",]);summary(m_c)
# Agree with algorithm
m_a_a <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
            data = tess2[tess2$treatment == "a_a",]);summary(m_a_a)
# Disagree with algorithm
m_a_d <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
            data = tess2[tess2$treatment == "a_d",]);summary(m_a_d)
# Agree with human
m_h_a <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
            data = tess2[tess2$treatment == "h_a",]);summary(m_h_a)
# Disagree with human
m_h_d <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
            data = tess2[tess2$treatment == "h_d",]);summary(m_h_d)
# Agree with both
m_ha_a <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
             data = tess2[tess2$treatment == "ha_a",]);summary(m_ha_a)
# Disagree with both
m_ha_d <- lm(jblame ~ female + age + pid7 + incomec + white + educ5, 
             data = tess2[tess2$treatment == "ha_d",]);summary(m_ha_d)

# Calculate ATMEs for Female
fem_a_a_coef <- coef(m_a_a)[2] - coef(m_c)[2]
fem_a_a_se <- sqrt(summary(m_a_a)$coefficients[2,2]^2 + summary(m_c)$coefficients[2,2]^2)

fem_a_d_coef <- coef(m_a_d)[2] - coef(m_c)[2]
fem_a_d_se <- sqrt(summary(m_a_d)$coefficients[2,2]^2 + summary(m_c)$coefficients[2,2]^2)

fem_h_a_coef <- coef(m_h_a)[2] - coef(m_c)[2]
fem_h_a_se <- sqrt(summary(m_h_a)$coefficients[2,2]^2 + summary(m_c)$coefficients[2,2]^2)

fem_h_d_coef <- coef(m_h_d)[2] - coef(m_c)[2]
fem_h_d_se <- sqrt(summary(m_h_d)$coefficients[2,2]^2 + summary(m_c)$coefficients[2,2]^2)

fem_ha_a_coef <- coef(m_ha_a)[2] - coef(m_c)[2]
fem_ha_a_se <- sqrt(summary(m_ha_a)$coefficients[2,2]^2 + summary(m_c)$coefficients[2,2]^2)

fem_ha_d_coef <- coef(m_ha_d)[2] - coef(m_c)[2]
fem_ha_d_se <- sqrt(summary(m_ha_d)$coefficients[2,2]^2 + summary(m_c)$coefficients[2,2]^2)

# Calculate ATMEs for Age
age_a_a_coef <- coef(m_a_a)[3] - coef(m_c)[3]
age_a_a_se <- sqrt(summary(m_a_a)$coefficients[3,2]^2 + summary(m_c)$coefficients[3,2]^2)

age_a_d_coef <- coef(m_a_d)[3] - coef(m_c)[3]
age_a_d_se <- sqrt(summary(m_a_d)$coefficients[3,2]^2 + summary(m_c)$coefficients[3,2]^2)

age_h_a_coef <- coef(m_h_a)[3] - coef(m_c)[3]
age_h_a_se <- sqrt(summary(m_h_a)$coefficients[3,2]^2 + summary(m_c)$coefficients[3,2]^2)

age_h_d_coef <- coef(m_h_d)[3] - coef(m_c)[3]
age_h_d_se <- sqrt(summary(m_h_d)$coefficients[3,2]^2 + summary(m_c)$coefficients[3,2]^2)

age_ha_a_coef <- coef(m_ha_a)[3] - coef(m_c)[3]
age_ha_a_se <- sqrt(summary(m_ha_a)$coefficients[3,2]^2 + summary(m_c)$coefficients[3,2]^2)

age_ha_d_coef <- coef(m_ha_d)[3] - coef(m_c)[3]
age_ha_d_se <- sqrt(summary(m_ha_d)$coefficients[3,2]^2 + summary(m_c)$coefficients[3,2]^2)

# Calculate ATMEs for pid7
pid_a_a_coef <- coef(m_a_a)[4] - coef(m_c)[4]
pid_a_a_se <- sqrt(summary(m_a_a)$coefficients[4,2]^2 + summary(m_c)$coefficients[4,2]^2)

pid_a_d_coef <- coef(m_a_d)[4] - coef(m_c)[4]
pid_a_d_se <- sqrt(summary(m_a_d)$coefficients[4,2]^2 + summary(m_c)$coefficients[4,2]^2)

pid_h_a_coef <- coef(m_h_a)[4] - coef(m_c)[4]
pid_h_a_se <- sqrt(summary(m_h_a)$coefficients[4,2]^2 + summary(m_c)$coefficients[4,2]^2)

pid_h_d_coef <- coef(m_h_d)[4] - coef(m_c)[4]
pid_h_d_se <- sqrt(summary(m_h_d)$coefficients[4,2]^2 + summary(m_c)$coefficients[4,2]^2)

pid_ha_a_coef <- coef(m_ha_a)[4] - coef(m_c)[4]
pid_ha_a_se <- sqrt(summary(m_ha_a)$coefficients[4,2]^2 + summary(m_c)$coefficients[4,2]^2)

pid_ha_d_coef <- coef(m_ha_d)[4] - coef(m_c)[4]
pid_ha_d_se <- sqrt(summary(m_ha_d)$coefficients[4,2]^2 + summary(m_c)$coefficients[4,2]^2)

# Calculate ATMEs for income
inc_a_a_coef <- coef(m_a_a)[5] - coef(m_c)[5]
inc_a_a_se <- sqrt(summary(m_a_a)$coefficients[5,2]^2 + summary(m_c)$coefficients[5,2]^2)

inc_a_d_coef <- coef(m_a_d)[5] - coef(m_c)[5]
inc_a_d_se <- sqrt(summary(m_a_d)$coefficients[5,2]^2 + summary(m_c)$coefficients[5,2]^2)

inc_h_a_coef <- coef(m_h_a)[5] - coef(m_c)[5]
inc_h_a_se <- sqrt(summary(m_h_a)$coefficients[5,2]^2 + summary(m_c)$coefficients[5,2]^2)

inc_h_d_coef <- coef(m_h_d)[5] - coef(m_c)[5]
inc_h_d_se <- sqrt(summary(m_h_d)$coefficients[5,2]^2 + summary(m_c)$coefficients[5,2]^2)

inc_ha_a_coef <- coef(m_ha_a)[5] - coef(m_c)[5]
inc_ha_a_se <- sqrt(summary(m_ha_a)$coefficients[5,2]^2 + summary(m_c)$coefficients[5,2]^2)

inc_ha_d_coef <- coef(m_ha_d)[5] - coef(m_c)[5]
inc_ha_d_se <- sqrt(summary(m_ha_d)$coefficients[5,2]^2 + summary(m_c)$coefficients[5,2]^2)

# Calculate ATMEs for white
white_a_a_coef <- coef(m_a_a)[6] - coef(m_c)[6]
white_a_a_se <- sqrt(summary(m_a_a)$coefficients[6,2]^2 + summary(m_c)$coefficients[6,2]^2)

white_a_d_coef <- coef(m_a_d)[6] - coef(m_c)[6]
white_a_d_se <- sqrt(summary(m_a_d)$coefficients[6,2]^2 + summary(m_c)$coefficients[6,2]^2)

white_h_a_coef <- coef(m_h_a)[6] - coef(m_c)[6]
white_h_a_se <- sqrt(summary(m_h_a)$coefficients[6,2]^2 + summary(m_c)$coefficients[6,2]^2)

white_h_d_coef <- coef(m_h_d)[6] - coef(m_c)[6]
white_h_d_se <- sqrt(summary(m_h_d)$coefficients[6,2]^2 + summary(m_c)$coefficients[6,2]^2)

white_ha_a_coef <- coef(m_ha_a)[6] - coef(m_c)[6]
white_ha_a_se <- sqrt(summary(m_ha_a)$coefficients[6,2]^2 + summary(m_c)$coefficients[6,2]^2)

white_ha_d_coef <- coef(m_ha_d)[6] - coef(m_c)[6]
white_ha_d_se <- sqrt(summary(m_ha_d)$coefficients[6,2]^2 + summary(m_c)$coefficients[6,2]^2)

# Calculate ATMEs for educ
educ_a_a_coef <- coef(m_a_a)[7] - coef(m_c)[7]
educ_a_a_se <- sqrt(summary(m_a_a)$coefficients[7,2]^2 + summary(m_c)$coefficients[7,2]^2)

educ_a_d_coef <- coef(m_a_d)[7] - coef(m_c)[7]
educ_a_d_se <- sqrt(summary(m_a_d)$coefficients[7,2]^2 + summary(m_c)$coefficients[7,2]^2)

educ_h_a_coef <- coef(m_h_a)[7] - coef(m_c)[7]
educ_h_a_se <- sqrt(summary(m_h_a)$coefficients[7,2]^2 + summary(m_c)$coefficients[7,2]^2)

educ_h_d_coef <- coef(m_h_d)[7] - coef(m_c)[7]
educ_h_d_se <- sqrt(summary(m_h_d)$coefficients[7,2]^2 + summary(m_c)$coefficients[7,2]^2)

educ_ha_a_coef <- coef(m_ha_a)[7] - coef(m_c)[7]
educ_ha_a_se <- sqrt(summary(m_ha_a)$coefficients[7,2]^2 + summary(m_c)$coefficients[7,2]^2)

educ_ha_d_coef <- coef(m_ha_d)[7] - coef(m_c)[7]
educ_ha_d_se <- sqrt(summary(m_ha_d)$coefficients[7,2]^2 + summary(m_c)$coefficients[7,2]^2)

# Create Tibble of results
mod_tib <- tibble(Treatment = rep(c("Agrees with algorithm", "Disagrees with algorithm",
                                    "Agrees with human", "Disagrees with human",
                                    "Agrees with algorithm and human", "Disagrees with algorithm and human"),7),
                  Moderator = c(rep("Trust in automation", 6), rep("Female", 6),
                                rep("Age", 6), rep("Ideology", 6),
                                rep("Income", 6), rep("White", 6),
                                rep("Education", 6)),
                  Effect = c(tia_a_a_coef, tia_a_d_coef, tia_h_a_coef, tia_h_d_coef, tia_ha_a_coef, tia_ha_d_coef,
                             fem_a_a_coef, fem_a_d_coef, fem_h_a_coef, fem_h_d_coef, fem_ha_a_coef, fem_ha_d_coef,
                             age_a_a_coef, age_a_d_coef, age_h_a_coef, age_h_d_coef, age_ha_a_coef, age_ha_d_coef,
                             pid_a_a_coef, pid_a_d_coef, pid_h_a_coef, pid_h_d_coef, pid_ha_a_coef, pid_ha_d_coef,
                             inc_a_a_coef, inc_a_d_coef, inc_h_a_coef, inc_h_d_coef, inc_ha_a_coef, inc_ha_d_coef,
                             white_a_a_coef, white_a_d_coef, white_h_a_coef, white_h_d_coef, white_ha_a_coef, white_ha_d_coef,
                             educ_a_a_coef, educ_a_d_coef, educ_h_a_coef, educ_h_d_coef, educ_ha_a_coef, educ_ha_d_coef), 
                  SE = c(tia_a_a_se, tia_a_d_se, tia_h_a_se, tia_h_d_se, tia_ha_a_se, tia_ha_d_se,
                         fem_a_a_se, fem_a_d_se, fem_h_a_se, fem_h_d_se, fem_ha_a_se, fem_ha_d_se,
                         age_a_a_se, age_a_d_se, age_h_a_se, age_h_d_se, age_ha_a_se, age_ha_d_se,
                         pid_a_a_se, pid_a_d_se, pid_h_a_se, pid_h_d_se, pid_ha_a_se, pid_ha_d_se,
                         inc_a_a_se, inc_a_d_se, inc_h_a_se, inc_h_d_se, inc_ha_a_se, inc_ha_d_se,
                         white_a_a_se, white_a_d_se, white_h_a_se, white_h_d_se, white_ha_a_se, white_ha_d_se,
                         educ_a_a_se, educ_a_d_se, educ_h_a_se, educ_h_d_se, educ_ha_a_se, educ_ha_d_se)) %>%
  mutate(EHI = Effect + 2*SE,
         ELO = Effect - 2*SE,
         Treatment = factor(Treatment, levels = c("Disagrees with algorithm and human",
                                                  "Agrees with algorithm and human",
                                                  "Disagrees with human",
                                                  "Agrees with human",
                                                  "Disagrees with algorithm",
                                                  "Agrees with algorithm")))
write_csv(mod_tib, here("Tables", "Study 1 moderation.csv"))


pt((-0.110/0.0574), df = 82, lower.tail = T) *2
pt((-0.0245/0.0564), df = 82, lower.tail = T) *2
pt((0.0381/0.0114), df = 342, lower.tail = F) *2
pt((0.0287/0.0116), df = 342, lower.tail = F) *2
pt((0.0304/0.0117), df = 342, lower.tail = F) *2

ATME <- ggplot(mod_tib) +
  geom_linerange(aes(x = Treatment, ymin = ELO, ymax = EHI, color = Treatment), size = 1) +
  geom_point(aes(x = Treatment, y = Effect), size = 2) +
  facet_wrap(~ Moderator, scales = "free_x") + coord_flip() +
  scale_color_npg() + theme_bw() + theme(legend.position = "none", axis.text = element_text(size = 8)) +
  labs(x = "Treatment", y = "ATME") +
  geom_hline(yintercept = 0)
ATME


####
# Add in Cohen's d for size of effects
####


temp_d <- tess2 %>%
  filter(treatment == "j" | treatment == "a_a") %>%
  mutate(treatment = factor(treatment, levels = c("a_a","j")))
cd_a_a <- cohen.d(temp_d$jblame, temp_d$treatment, na.rm = T)
cd_a_a
cd1 <- c("Alg agree", cd_a_a$estimate, cd_a_a$conf.int)

temp_d <- tess2 %>%
  filter(treatment == "j" | treatment == "a_d") %>%
  mutate(treatment = factor(treatment, levels = c("a_d","j")))
cd_a_d <- cohen.d(temp_d$jblame, temp_d$treatment, na.rm = T)
cd_a_d
cd2 <- c("Alg disagree", cd_a_d$estimate, cd_a_d$conf.int)

temp_d <- tess2 %>%
  filter(treatment == "j" | treatment == "h_a") %>%
  mutate(treatment = factor(treatment, levels = c("h_a","j")))
cd_h_a <- cohen.d(temp_d$jblame, temp_d$treatment, na.rm = T)
cd_h_a
cd3 <- c("Hum agree", cd_h_a$estimate, cd_h_a$conf.int)

temp_d <- tess2 %>%
  filter(treatment == "j" | treatment == "h_d") %>%
  mutate(treatment = factor(treatment, levels = c("h_d","j")))
cd_h_d <- cohen.d(temp_d$jblame, temp_d$treatment, na.rm = T)
cd_h_d
cd4 <- c("Hum disagree", cd_h_d$estimate, cd_h_d$conf.int)

temp_d <- tess2 %>%
  filter(treatment == "j" | treatment == "ha_a") %>%
  mutate(treatment = factor(treatment, levels = c("ha_a","j")))
cd_ha_a <- cohen.d(temp_d$jblame, temp_d$treatment, na.rm = T)
cd_ha_a
cd5 <- c("Alg & hum agree", cd_ha_a$estimate, cd_ha_a$conf.int)

temp_d <- tess2 %>%
  filter(treatment == "j" | treatment == "ha_d") %>%
  mutate(treatment = factor(treatment, levels = c("ha_d","j")))
cd_ha_d <- cohen.d(temp_d$jblame, temp_d$treatment, na.rm = T)
cd_ha_d
cd6 <- c("Alg & hum disagree", cd_ha_d$estimate, cd_ha_d$conf.int)

cd_tot <- data.frame(rbind(cd1,cd2,cd3,cd4,cd5,cd6)) %>%
  rename(treatment = V1, cd = V2, cdl = lower, cdh = upper) %>%
  mutate(treatment = factor(treatment, levels = c("Alg & hum disagree",
                                                  "Alg & hum agree",
                                                  "Hum disagree",
                                                  "Hum agree",
                                                  "Alg disagree",
                                                  "Alg agree")),
         cd = as.numeric(cd),
         cdl = as.numeric(cdl),
         cdh = as.numeric(cdh)) %>%
  ggplot() +
  geom_bar(aes(x = treatment, y = cd), stat = "identity") +
  geom_errorbar(aes(x = treatment, ymin = cdl, ymax = cdh), width = 0.4, alpha = 0.9, size = 1.3) +
  coord_flip() +
  labs(y = "Cohen's d", x = "Treatment") +
  scale_y_continuous(limits = c(0, 0.5, 0.1)) +
  theme_bw()
